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ABSTRACT 

An interactive computer package for generating a finite 
element mesh has been developed as a preprocessor for two 
and three-dimensional aquifer models. This scheme has the 
capability to generate triangular and quadrilateral elements 
in two-dimensions and tetrahedral and parallelepiped 
elements in three-dimensions. It is based on the 
Isoparametric mapping of quadrilaterals and can be used on 
flow domains that are described by both regular and 
irregular boundaries. The system is completely user 
interactive, utilizing a digitizing tablet and a colour 
graphics terminal to involve the user in all phases of the 
mesh generation, and also decrease the time and effort that 
is involved in more traditional methods of mesh 
construction. There are facilities to incorporate 
pre-positioned pumping wells and observation wells into the 
mesh, aS well as constant head nodes, boundary fluxes and 
leakance elements. Output from the package is stored in a 
file in a format that is easily adaptable to most popular 
finite element codes. Some features of the code include: the 
ability to produce a hard copy of the mesh in any size 
within the limitations of the plotting hardware, and an 
optimization algorithm that renumbers nodes so as to 
increase computational efficiency. In addition to a 
description of the code, guidelines and hints are given for 
dealing with unusual shaped domains and for improving mesh 


design, and examples of the application of the package to a 
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1. INTRODUCTION 

Finite element methods for solving groundwater flow 
problems require the discretization of the flow domain, 
which involves the partitioning of the system under 
consideration into a grid or a mesh. Typically the mesh is 
constructed by hand through drawing a grid over the flow 
domain and manually entering the mesh parameters into a data 
set. This procedure, especially when one encounters more 
complex flow systems is very tedious, time consuming and 
Subject to errors (Klienstreurer, 1980). 

The need to automate this process is therefore apparent 
(Narasimhan and Witherspoon, 1980, and Klienstreurer, 1980). 
Techniques based on point digitizers have proven useful in 
reducing the time for data entry of a manually prepared 
grid. However, there has been relatively little work in 
adapting fully automatic methods used for example in 
Structural engineering to hydrogeologic problems. These 
schemes have been widely used, because they provide the 
ability to quickly and accurately change the density of the 
mesh to provide an optimal solution to a given problem. 

In the fields of structural and automotive engineering 
the problem of automatic mesh generation appears somewhat 
less complex than for hydrogeological applications because 
of regular external and internal geometries. In 
hydrogeologic applications, mesh generation is complicated 
by the irregular geometry of the flow domain. Specifically 
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-complex internal geometries, where it is necessary to 
model the boundaries of lithologic units; 
-complex external geometries, as for example in the case 
of buried channels; 
-providing specific locations for observation and 
pumping wells that must correspond to nodal locations; 
and 
-the necessity of compressing the grid around a pumping 
well to increase the mathematical accuracy in an area 
where the hydraulic head may be changing quite markedly. 
Such constraints on the mesh to be generated are 
usually not considered in the schemes presently available 
today. Existing techniques applied to groundwater modelling 
have been largely problem or program specific 
(Klienstreurer, 1980). These schemes generally involve 
fitting the flow domain to an idealized grid of regularly 
positioned elements. This kind of constraint usually makes 
it difficult to accurately model an irregular flow domain. 
A more general procedure for mesh generation ina 
pre-processor form has been developed by Sartori et al. 
(1982). The meshes illustrated in their work appear to be 
accurate and well conditioned. This scheme fits a regularly 
positioned grid upon an irregular geometry such that 
elemental boundaries do not necessarily correspond to 
transmissivity boundaries. 
Thus, there is still a need in hydrogeological 


applications to extend existing procedures for mesh 
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generation. The proven techniques in the field of structural 

engineering provide the logical starting point. Although 

many schemes exist for mesh generation, they are generally 
based on the following basic philosophies: 

1. the mesh pattern is established by the computer from a 
minimal amount of information supplied in digital form; 
and 

2. the positioning of the mesh is established by 
graphical-computer interaction using plotting tables or 
direct light pen displays (Zienkiewicz and Phillips, 
1971). 

The goal of my work is to extend these basic 
philosophies of mesh generation to hydrogeological problems. 
Specific objectives are: (1) to develop a mesh generator, 
capable of constructing a network of elements in two- or 
three-dimensions a level of complexity appropriate for most 
problems encountered in aquifer hydraulics, (2) to 
demonstrate the usefulness of this approach in solving 
practical problems, and (3) to demonstrate how grid design 
can improve finite element simulations. 

The scheme for mesh generation used in this work is 
discussed in Chapter 2. It is based on an Isoparametric 
method developed by Zienkiewicz and Phillips (1971). The 
System is designed for maximum flexibility, user 
interaction, and is capable of generating triangular and 
quadrilateral shaped elements in a plane and projecting 


these elements in the third dimension to produce tetrahedral 
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and parallelepiped (brick type) elements. The mesh can be 
layered and hence can account for the patterns of geologic 
layering with relative ease. 

The demonstration of the mesh generation scheme 
provided in Chapter 2 involves the simulation of hydraulic 
head distributions for both a two- and three-dimensional 
problem. The mesh for the two-dimensional case is prepared 
both manually and with the mesh generation program. The 
results of these two Simulations are compared to verify the 
automatic procedure and to demonstrate its utility. The 
three-dimensional simulation is a more complex case designed 
to show how the mesh generation procedure can be applied to 


layered geological systems. 
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2. DESIGN AND IMPLEMENTATION OF A FINITE ELEMENT MESH 


GENERATOR 


2.1 EXISTING METHODS 

_ Techniques for automatic mesh generation have been in 
existence in their most basic form since the late 1960's. 
Thwacker (1980) provides an exhaustive review of these mesh 
generation techniques, including an outline of the basic 
theory, and an analysis of the drawbacks of each scheme. 
Current mesh generators employ one or a combination of two 
of the following schemes: 
te automatic triangulation; 
2. coordinate transformations; 
3. smoothing procedures; or 
4, blending functions (Yerry and Shephard, 1983). 

The automatic triangulation technique, a scheme 
demonstrated in the works of Cavendish (1974), Sadeck 
(1980), and Suhara and Fukuda (1972), allows the greatest 
flexibility in modelling a very general region. Cavendish 
(1974) uses a random number generator to fill the 
computational domain with nodes to some specified density. 
From this population of nodes, nodes are selected to form 
elements in a well conditioned mesh. Nodes that are not 
selected are extraneous and ignored. 

Sadeck (1980) takes this approach a step further. 
Instead of randomly generating the coordinates of the nodes, 


the nodes are located on the basis of the position of 
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regularly spaced nodes around the boundary of the region in 
which the mesh is being constructed to form a well 
conditioned mesh. This procedure eliminates the need to 
generate the nodes first and then sort through them to 
produce the best conditioned element from the distribution 
of nodes. 

The method used by Sadeck (1980) works well for regular 
geometric shapes, or simple forms with irregular geometry. 
For regions with complex geometry, which are discretized 
uSing quadrilateral elements, this method tends to form 
poorly conditioned elements especially where the boundary 
forms an acute angle. My experience with this method has 
shown that, although relatively simple to program it does 
require a substantial number of support subroutines to 
handle complex shapes. Yerry and Shephard (1983) state that 
although the triangulation algorithms have been developed 
that will work for general two-dimensional shapes, they have 
not proven popular because they occasionally produce 
ill-conditioned meshes and in many cases require substantial 
user interaction. 

Coordinate transformations are used to map 
geometrically complex figures to a regular grid in the 
transformed domain. Having defined the transformation, the 
inverse transformation can be used to revert to the original 
domain. A disadvantage of this technique is that although 
the mesh is apparently well conditioned, depending on the 


coordinate transformation the mesh in reality may be very 


of Bean sf essenimiis eve Rees aie. ee 


o¢ Mets dpuowts sice neds brs sani? 298er sad ® 

noi sudigse lh SAw mot? soemsls Sstois { baes sasd ody 
. F : : 

siuges 103 Lie erro" (O82: 1 79862 qd bau howsom 

-vasoneae t6fupets dst 4 29702 -@pgare ts ,2eqadet 

bax iteis2i6 976 perry yrismeag <signe> atin ers. 

org? o? ebhns: bobs pili: alagnte fay fe vetelivieep pnd 


yisbrucd ody sicitw yilsissces eteamnsis beaters — —_ 


fare siete (SBE) Brediqaed2 bre yu Tsay Bogan akauabiiets 


le ' E 


beqoieveb need ovat anivizogta rotisivenséas ony squety 


ved yous ,eeqgete Lanciengmib-or2 isisase 10 ; Atow itiwg 

soubeid. yiienoienzz0 ysis 4 seusoss jaLugag navorg 4 

ey. A : ; > 

giinscedve etivpss 2es26> yen ne bas soneea etia’ se 
| = 


or : ». 3 hs ’ 2 A "i igiseereant yy sf 


- * : 


qem og peev ate eno? saute eaee! 
en nt bite. retupby a ot: easupit ete 
nis 108 suereenian on? oo he — 
saigize wf et ions oF. heme ae meet ver 


os 


pode “= f 


fag) 


distorted. As well, another disadvantage of this scheme is 
that the coordinate transformation must be provided before 
the computation can begin. This requirement is analagous to 
constructing an irregular grid (Thwacker,1980). 
Isoparametric methods (Zienkiewicz and Phillips, 1971) 
or, smoothing procedures as Yerry and Shephard (1983) 
describe them, involve the smooth distribution of nodes 
within the computational domain. The positioning of the 
interior nodes is determined by the mean of the coordinates 
of 1ts neighbours. This condition is equivalent to an 
optimal grid, which corresponds to the equilibrium condition 
of a system of springs connecting the nodal points 
(Thwacker, 1980). This scheme requires that opposite 
boundaries of the domain have the same number. of elements. 
This constraint is a disadvantage because the user can not 


easily change the density of the elements within the domain. 
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2.2 DESIGN OF A FINITE ELEMENT MESH GENERATOR FOR TWO- AND 
THREE-DIMENSIONAL PROBLEMS 

My design for an automatic mesh generator for 
three-dimensional groundwater flow problems is based upon 
the work of Zienkiewicz and Phillips (1971), which involves 
the application of isoparametric coordinates, or smoothing 
procedures as Yerry and Shephard (1983) had described them 
in the previous section. The work by Zienkiewicz and 
Phillips (1971) deals mainly with two-dimensional meshes, 
and indicates that the method can be easily adapted to 
three-dimensional problems. 

The code that I have developed has the capability of 
generating finite element meshes in both two- and 
three-dimensions. It discretizes a domain using either 
quadrilateral or triangular elements in two-dimensions, and 
parallelepiped or tetrahedral elements in three-dimensions. 

To facilitate both the rapid and simple entry of data 
into the code, a menu oriented approach has been implemented 
using a Summagraphics Corp. ID series tablet digitizer to 
input the bulk of the digital information. This approach 
greatly simplifies digital input of coordinate information, 
and provides the user with flexibility in controlling the 
pattern of mesh generation. 

Coupled to the digitizer is an Advanced Electronics 
Design AED512 colour graphics terminal that completely 
Supports the graphics packages necessary to display and 


manipulate the networks that are generated. The Integrated 
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Graphics (Davis, 1977) system is a library of 
FORTRAN-callable subroutines for performing graphics 
operations. These systems are all supported by the Amdahl 
5860 mainframe at the University of Alberta. 

With this hardware, the amount of data that must be 
input manually can be greatly minimized to generate a mesh 
consistent with the following criteria: 

1. the geometry of the boundaries, whether singly or 
multiply connected, must be accurately described; 

2. a facility for describing the spatial variability in the 
hydraulic properties of the geologic material, wherever 
these occur; 

3. a facility for grading the mesh so as to achieve 
different kinds of mesh compression in specific 
Sub-areas of the domain; 

4. an achievement of shapes which are approximately 
equidimensional and do not lead to numerical 
pll=conditiondangs 

5. a renumbering system which results in the best 
computational efficiency (Zienkiewicz and Phillips, 
POM Co) Pins. 1973): 

The boundary of the domain or sub-domain in any plan 
view can be adequately described by digitizing points around 
this boundary. Zones of differing transmissivities are 
digitized separately, with the program establishing the 
required continuity in the mesh geometry along boundaries. 


The manner in which the other three criteria are implemented 
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in the procedure for mesh generation will become apparent in 


the following sections. 


The mesh generating scheme is based on the 


isoparametric mapping of quadrilaterals. For any shaped 


domain the coordinates of the interior nodes can be 


described by the following equations (Zienkiewicz and 


Phillips, 1971): 


(1) 


the shape factor determining the placement 


the nodes around the domain boundary. 


X=N ; Xj; 
yY=Niyi 
Z=N ,Z; 
where: N;, = 
X; = x-coordinate 
being mapped 
Yi = y-coordinate 
being mapped 
Z, = zwzrcoordinate 


being mapped 


ofra 
onto 


of a 
onto 


ofaa 
onto 


regularly spaced node 
the computational domain. 


regularly spaced node 
the computational domain. 


regularly spaced node 
the computational domain. 


The shape factor N, is described in terms of a 


coordinate transformation based upon the geometry of the 


domain boundary. In the field of hydrogeology, the geometry 


of the flow domain is not easily described as a simple 


equation. 
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The approach taken in this thesis involves defining the 
boundaries of the figure by a series of regularly spaced 
nodes. The spacing of these nodes is based on the desired 
distribution of elements within the figure. In 
three-dimensions the nodes are regularly spaced along the 
twelve edges or folds of the figure. 

To obtain a regular spacing along an irregular 
boundary, a method based on Lagrange polynomials (Hornbeck, 
1975) is incorporated into the code. Essentially this 
Lagrangian technique produces polynomials that are "tied 
down" (Hornbeck 1975) at the data points given in the 
digitization and interpolates positions in between. In the 
code, boundaries of the domain or sub-domain are represented 
by equations of up to a third order polynomial, which can 
represent most geometries that may be encountered. 

The Lagrangian method interpolates evenly spaced points 
from an irregular distribution of sample points. For this 
method, consider a series of points described by a function 
f(x,;) where the x; are in general not evenly spaced and i 
can take on integer values 0 to n (meaning there are nt+1 
data points). To fit a polynomial interpolation between the 
x,;'S 1S to fit a polynomial of degree n to these nt+1 points. 
The mechanics of this method in this code is a moving 
average approach through the data points as a constantly 
changing polynomial. That is, utilizing four data points at 
a time, the code calculates evenly spaced points within a 


range represented by these points. As the point generation 
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progresses, the sample points change to keep the generated 
points within the range represented represented by the four 
data points. 

Because this scheme is based on the mapping of 
quadrilaterals, it is necessary to designate four corner 
points (in three-dimensions, eight corner points) of the 
irregular figure to obtain an equal number of regularly 
spaced nodes on opposite sides of the figure. With the 
digitized figure displayed on the graphics terminal, the 
user, by positioning the cursor defines the four (or eight) 
corners of the figure. In dealing with irregular geometries, 
corner points are not immediately obvious but must be 
designated in a manner to provide an even distribution of 
elements as illustrated in Figure 1. 

With the boundary nodes evenly spaced around the 
boundary of the domain under consideration, the boundaries 
are completely defined but not in the equation form as in 
the work of Zienkiewicz and Phillips (1971). The 
distribution of the interior nodes becomes based on the 
distribution of the boundary nodes and may be described by 


the following set of equations: 


N,=(x,-x2)/NELM, 
Ny=(yi-y2)/NELM, 


N,=(2,-Z2)/NELM, 23 
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where: N, = the factor determining the distribution 
of the nodes in the x-direction. 


N, = the factor determining the distribution 
of the nodes in the y-direction. 


N, = the factor determining the distribution 
Sieune nodes jimtbesz-direction. 


xX, & X2 = corresponding x-coordinates of nodes 
on opposite sides of the domain. 


yi & X2 = corresponding y-coordinates of nodes 
on opposite sides of the domain. 


z; & Z2 = corresponding z-coordinates of nodes 
on opposite sides of the domain. 


NELM, = the number of elements to be generated 
in the x-direction. 


NELM, = the number of elements to be generated 
generated in the y-direction. 


NELM, = the number of elements to be generated 
in the z-direction. 


As the nodes are’ being generated, the elemental incidences 


are being described simultaneously. 


2.3 DESCRIPTION OF THE PROGRAM 

The mesh generation program consists of 13 parts: the 
main program and 12 subroutines: SMOOTH, DSTRB, NODE, ELMNT, 
SETUP, OPTNUM, NDEGEN, INTANG, ASKPLT, CNCTQ, CNCTT1, 
CNCTT2, as well as the Integrated Graphics Library 
Subroutines. The program is presently dimensioned to 
generate 5000 nodes and 3000 elements, which should be more 


than sufficient for most problems. 
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MENU 


Generate a Mesh(after digitizing) 
Output Data for Later Use 

‘Input Boundary Fluxes 

Future Use 

Digitize a Transmissivity Boundary 
Digitize an Observation Well 
iotieriieae a Pumping Well 

“Pic’ Constant Head Nodes 

Pic’ Leakance element 

Remove Last Transmissivity Zone Generated 
Remove Last Point Digitized 
Stop/Begin Again 


Re-display Menu 


Figure 3. Menu illustrating options available to a user to generate a mesh. 
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The MAIN program initializes all the indices used in 
the entire program, reads in the title of the simulation and 
interactively sets the baSic options of the program. These 
options determine: whether one is working in English or 
metric units; the scale; whether it is to be a 
two-dimensional or three-dimensional simulation; type of 
elements required; and if in three-dimensions whether a 
Steady state or transient simulation is to be performed. The 
MAIN program also contains the Menu from which any of the 
other subroutines will be called. This is implemented as a 
Series of 13 switches that may be accessed at any time 
during the construction of the mesh in a manner outlined in 
the flow chart, Figure 4. The 13 switches allow the user to 

perform a variety of tasks to synthesize a mesh. 

Switch No. 1 simply re-displays the menu. This feature 
was incorporated into the design of the program to provide 
the user with an easy way to recall all of the options. 
Carrying most of the operating instructions with the program 
eliminates the need for a paper copy of the menu that could 
be misplaced. 

Switch No. 2 activates a built-in procedure to restart 
Or terminate the program. The instructions outlined in the 
user's guide in Appendix A explains how to start the 
program. However because the user has options throughout the 
entire program that may never be used, an on/off switch is 


provided to terminate the program. 
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Switch No. 3 deletes the previously digitized point. 
This feature is incorporated so that one could take care of 
such hardware problems such as a faulty digitizer switch, 
Static on the line and human error. 

Switch No. 4 allows one to remove the mesh generated 
for the last transmissivity zone incorporated into the main 
mesh. This feature avoids repeating large amounts of data 
entry in the case where a mistake in the digitizing process 
goes unnoticed. 

Switch No. 5 allows the user to define leakance 
elements, (ie: recharge associated with a river). This 
feature only works in two-dimensions, uSing a graphic 
routine to pick out particular elements ("Pic" routine) with 
the. program prompting a numerical value of leakance. 

Switch No. 6 defines constant head nodes in 
two-dimensions through a "Pic" routine and in 
three-dimensions constant head boundary faces. The user 
cannot define individual constant head nodes in 
three-dimensions with this program because the "Pic" routine 
cannot easily define the z-coordinate. 

Switch No. 7 defines the location of a pumping well in 
the x-y plane. As well the program prompts the user to input 
information relating to the pump schedule. The position of 
the well is incorporated into the mesh as a node and in 
three-dimensions is continuous throughout all layers 


subsequent to its initial nodal definition. 
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Switch No. 8 defines the location of an observation 
well in a Similiar manner to the pumping well. 

Switch No. 9 allows the user to define the boundaries 
of a transmissivity zone with the program prompting for 
values of storativity and transmissivity (in 
three-dimensions, elemental hydraulic conductivities in x- 
,y-~ and z-directions). 

Switch No. 10 is designated for future use. It was 
Originally used to renumber the nodes to reduce bandwidth of 
the mesh, but this feature was subsequently incorporated 
into Switch No. 12. 

Switch No. 11 designates boundary fluxes across the 
entire face of a mesh. It was felt that a boundary flux 
paar was needed in the mesh generation scheme, but to 
"Pic”™ out individual elements would be too costly and time 
consuming. 

Switch No. i2 renumbers the nodes of the mesh to reduce 
bandwidth (Collins, 1973) and outputs the information in a 
format that is compatible with the Verge-Frind (Verge, 1972) 
two-dimensional code dealing with transient simulations for 
confined aquifers. In three-dimensions, a format compatible 
with Verge's (1975) three-dimensional groundwater flow model 
with minor changes in the input. The numerical information 
is output to device 7 and the plotting information to device 
9. The workings of Collins' (1973) renumbering scheme will 


be discussed in Section 4. 
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Switch No. 13 activates the mesh generation procedure. 
This procedure is interactive with the user, allowing him to 
define the configuration of the mesh and the elemental 
density. If any pumping wells are present, the program will 
concentrate elements around the well and smooth the 
distribution before it is displayed on the terminal. Having 
displayed the mesh on the terminal, the user can rotate it 
to better visualize its faces and plot it on a hardcopy 
plotter. 

The individual subroutines used in the code correspond 
to no particular switch discussed previously. A switch may 
have aS many as five subroutines age Fee Wilstiwel't sO7 dt 
may have none. For this reason, the function of each 
Subroutine will be discussed briefly. 

Subroutine SMOOTH, used in conjunction with subroutines 
CNCTQ, or CNCTT1, or CNCTT2, concentrates elements around 
pumping wells and smooths the distribution of elements in 
the vicinity of the wells. Subroutine CNCTQ applies to 
quadrilateral and "brick type" elements around wells, 
whereas CNCTT1 and CNCTT2 apply to triangular and 
tetrahedral elements. 

The concentrating subroutines CNCTQ, CNCTT1, or CNCTT2 
partition the previously generated elements around a pumping 
well into smaller elements. The smaller elements are 
mumbered and the incidences defined. The incidences of the 
elements around these partitioned elements also have their 


incidences redefined. A large part of these subroutines can 
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be thought of as being largely devoted to bookeeping duties 
in the storage of elemental incidences. 

Subroutine SMOOTH takes the coordinates of the newly 
generated nodes in CNCTQ, CNCTT1, or CNCTT2 and repositions 
them uSing an averaging technique based on the position of 
the surrounding nodes. This provides for a smooth 
distribution of elements in the vicinity of the pumping 
well. 

Subroutine DSTRB re-distributes the digitized points 
around a perimeter of the transmissivity zone using the 
Lagrange interpolation routine (Hornbeck, 1975) 

Subroutine NODE uses the redistributed nodes' positions 
to produce nodes within the interior of the transmissivity 
zone, with the elemental incidences being described by calls 
to subroutine ELMNT. 

Subroutine SETUP and OPTNUM are based on algorithms 
developed by Collins (1973) and modified for use in this 
program. SETUP, as the name implies, sets up arrays that 
describe the interconnection relationship of each node. The 
purpose of these arrays is to make it easier for the 
renumbering routine to know which nodes are interconnected, 
hence the differences in numbering these nodes should be 
small. OPTNUM takes the results from SETUP and renumbers the 
nodes in the mesh based on each node's connection with its 
neighbours so that the bandwidth is reduced. 

Subroutines INTANG and NDEGEN are used in conjunction 


with the subroutines for mesh compression. INTANG calculates 
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the interior angle formed by the positioning of three 
specified nodes. NDEGEN uses the angle determined in INTANG 
and a specified radial distance defines a node in 
rectangular coordinates from polar coordinates. 

Subroutine ASKPLT provides the user with the capability 
to generate a hardcopy plot. ASKPLT also allows the user to 
change the size of the plot. The default size is the same 


Size that is viewed on the screen. 
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3. APPLICATIONS OF THE FINITE ELEMENT MESH GENERATION 


PROGRAM 


3.1 TWO-DIMENSIONAL MESH GENERATION 

To illustrate the operation and ease of use of this 
program, the following sample problem is included. Figure 2 
depicts the two-dimensional domain of interest. The 
groundwater system consists of a confined aquifer enclosed 
on all sides by no flow boundaries, with a river providing a 
source for induced infiltration. There is one pumping well 
and three observation wells as shown on Figure 2. The 
transmissivity of the aquifer east of the river is 4.0x10~* 
mnevsyewhide wests offithebriversite ise 20x 0 #em7/siorhe 
pumping well is assumed to be pumped at a constant rate of 
6.32x10°° m?/s. The leakance rate from the river is 1.0x1078 
m/s. This is the basic information that is required to 
produce a mesh. 

Once execution of the program has begun using the run 
command in Appendix A, the first prompt will be to enter: 1) 
a title; 2) the units one is working in, (ie: English or 
metric); 3) the scale of the figure from which one will be 
producing the mesh; 4) whether one is working in two- or 
three-dimensions; 5) and whether one desires quadrilateral 
(brick type) or triangular (tetrahedral) type elements. 

Next, the user is requested to clear the digitizer, and 
digitize the lower left, lower right, and upper left corners 


of the page of the figure that is to be digitized. The 


24 


eine Yo wen. Zo sans bos sotdaneae. ‘ego stortautifd 


_bebutont et meitorq signee poiwel to? ait 
edt . deeaeant to A hina 5. nena a? 
besoLons aoreuipe beni iao>- 5: ¢d. ada fens medaye 9 
£ entbivexd 79yi2 6 024 cot sabntod: wot on yd ane 

noderehs iinet! beoubat ae 


eles noitevisede -_ 


sit <8 grup it 10 wore 2a 
*-Ofno.® gi tevir ode ro sens /yateags gid) Foo gs vee oy 
. opt S ith: tov ie “aria” to 7288 ati 4 aN ae ; 
ed Gt pemusae et 1 Laie c 


¢ sxupit 


iiew phigmugq sre 21 S198 


sat. 2 
vo “916% yemsea0° & 76 DecmUC 
Sivies ad mo s1573- scngkee | ot. eum “tors 
{emai nF sae anit. et war 


doom & 


o- “OF x0 . } 21 I 
. 63-5 +ivpet a tds 3. NOt 


igs ot eniay nuged / ae aex803g ait 20 aktpoene 9 
} tera OF: sditiw aqnio%e t@aiz’ heal pecan 
s¢ deriena dahl st eqtaow ak ‘and See a 

ad ifiw amo detdw nox. swe sad: 36 ‘oleae 


_ 0 Eig 


ae | 


F. 


20 


purpose of this step is to scale the figure for display on 
the graphics terminal. 

The menu is now displayed on the screen and the user 
follows the steps outlined in the flow chart, Figure 4. 

The first step is to digitize the boundaries of a given 
transmissivity zone, and fixing the transmissivity value of 
all the nodes in this zone. Once the transmissivity boundary 
is defined, the user would then re-display the menu by 
digitizing Switch No. 1. Before moving to the next 
transmisSivity zone, it is necessary to define all 
observation and pumping wells. For the problem depicted in 
Figure 2, the first transmissivity zone contains two 
observation wells and one pumping well. By selecting Switch 
No. 7, the user can enter the pumping information (ie: start 
and finish times of the simulation, number of time steps and 
pumping rates) as requested, and digitize the pumping well 
location. Observation well locations are input using Switch 
No. 8. With this information, the user is able to generate a 
mesh for this transmissivity zone by selecting Switch No. 
13. The outline of the transmissivity zone is displayed onto 
the screen. Because this scheme is based on the mapping of 
quadrilaterals, the user is prompted to indicate with the 
cursor the four corners of the figure (even if one desires 
triangular elements). 

Once the four corners are indicated, the program 
generates boundary nodes using a method of interpolation 


through Lagrange polynomials (Hornbeck, 1975). The prompts 
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request the four corners of the figure, (see Figure 1), and 
request the number of elements in both the x and y 
directions. 

An important feature of this program is that 
observation and pumping wells are digitized prior to the 
mesh generation so that nodes can be placed at these 
locations and incorporated in the mesh. 

If the transmissivity zone contains any pumping wells, 
the mesh is smoothed. This smoothing routine adjusts the 
positions of each of the interior nodes to make sure that 
each of the elements are well conditioned. Also, elements 
are concentrated around pumping wells (see Section 4) to 
mathematically handle the large changes in hydraulic head 
encountered here. The resultant mesh is then displayed a 
the screen with the query "PLOT?". Because we have not 
completed our mesh generation process the response would be 
,"no", and the menu would then be redisplayed on the screen. 

The next area of interest would be the river because 
the boundaries of elements must correspond to the boundaries 
of the river for the purposes of defining leakance elements. 
Selecting Switch No. 9 again, we enter the boundaries of the 
river aS a separate transmissivity zone, and generate a mesh 
as before. Leakance elements are defined by selecting Switch 
No. 5, which the user does by indicating, with the cursor, 
the leakance elements and the value of leakance. This 
routine can be rather tedious, but does allow one to define 


leakance elements anywhere on the mesh. 
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The user could continue the mesh generation process 
with the area west of the river. 

With the mesh completely drawn, the user can choose to 
have it plotted by sending the plot data into a file on 
Output device 9 for the hard copy plotter. 

By selecting Switch No. 12, the nodal points are 
renumbered to reduce bandwidth (Collins, 1973), and output 
into a file attached to device 7 in a format that is 
compatible with the Verge-Frind (Verge, 1972) 
two-dimensional, finite element program for confined aquifer 
Simulations. 

The mesh generated using this procedure is shown in 
Figure 5. It can be compared with the hand-drawn version 
illustrated in Figure 6. The hand-drawn mesh has 104 nodes, 
177 elements, with a bandwidth of 23, whereas the computer 
generated mesh has 160 nodes, 272 elements with a bandwidth 
of 21. The hand drawn version required approximately 6 hours 
to draw, number, enter the data file and check for errors. 
The computer generated mesh required approximately 6 seconds 
of processing time on an Amdahl 5860 mainframe to produce, 
and approximately 20 minutes at the terminal to enter the 
necessary data. 

The simulated hydraulic head distribution due to seven 
days of pumping, based on the computer generated mesh, is 
presented in Appendix C. By comparing these results with 
those simulated using the hand-generated mesh, it is 


possible to see that the results are very similiar. The 
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Original purpose of this excercise was to determine whether 
it would be worthwhile to go to the expense of installing 
observation well number 3. According to the results of the 
Simulations we can See that in both instances one would have 
to question the expense, for the drawdown indicated at this 
point would only be 0.03 m for the hand-drawn version, with 
a Similiar value of 0.02 m for the computer generated 
version. The values of the hydraulic head differ from less 
than 0.01 ft in the area of observation well 3 to 3 ft. 
around the pumping well. The discrepancies between the two 
Simulations are approximately in the order of 10%. The 
source of this discrepancy is not immediately apparent, but 
may be due to several sources. The sources could be 
accounted for in differences in mesh density, computational 
efficiency, and inherent inaccuracies within the simulation 
algorithm. 

Variations in mesh density could account for a certain 
amount of this discrepancy. The higher density may have been 
better able to cover the changes in flux. This feature will 
be more fully discussed in the following chapter. As well, 
there may be inherent problems within the simulation program 
itself. In my experience, as much variation as is seen 
between these two simulations can be seen when viewing many 
hand drawn versions. 

To further illustrate the capabilities of my algorithm, 
additional two-dimensional examples are included in this 


section. Figures 11a, 11b and 11c illustrate how easy it is, 
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with an automatic mesh generator, to increase the nodal 
density within a single rectangular domain. The facility to 
alter the density of the nodes within a domain is 
represented by a prompt within the mesh generation routine 
accessed through Switch No. 13. This prompt requests that 
the user indicate the number of elements desired in the x 
and y directions and, if working in three-dimensions, the 
z-direction. It is therefore the duty of the user to insure 
that the resultant elements are as well conditioned as 
possible. 

Figures 8 through 14 are included to show what meshes 
actually look like when the domains take on a somewhat more 
complex shape. Figure 13 represents a shape that may be 
encountered if one were to model a buried channel situation. 
The four corners were placed as shown with an "X", Figures 
13 and 14 show how important the initial placement of four 
corners is in controlling the character of the grid. With a 
triangular shaped region, it can be seen that the mesh in 
Figure 14 represents a smooth distribution of elements from 
the apex to the opposite side. In Figure 14, with the corner 
nodes placed as shown, there is a stretching of the 
elements, producing some ill-conditioned elements at the 
lower boundary. This last figure illustrates how it is 
possible to alter the mesh density along a given boundary by 
altering what would be an obvious placement of corner nodes 
to one that would alter the distribution of nodes along a 


given boundary. In this way any adjacent zones would not be 
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constrained to a set mesh density. 


3.2 THREE-DIMENSIONAL MESH GENERATION 

To illustrate the capabilities of the code in 
three-dimensional applications, consider the following 
hypothetical problem of assessing the influence of an open 
pit mine on groundwater levels. The map illustrating the 
layout of the simulation is shown in Figure 7. For this 
problem the system is discretized into eight transmissivity 
zones representing a layered sedimentary sequence striking 
approximately east-west and dipping southward. Layers E,G 
and H are assumed to outcrop at the surface. The 


hydrogeological parameters are outlined in Table 1. 


TABLE 1. 


SUMMARY OF BASIC DATA FOR INDIVIDUAL LAYERS IN THE MODEL 


Layer Predominant Thickness Hydraulic Trans. 
No. Lithology (m) Cond. (m/s) (m?/s) 

H limestone 125 Tx 10™* 1.25x@ 07 * 

G limestone 30 Pxi1077 30x07 ° 

F shale 20 tx 05s 2.0xH07? 

E shale 15 ct 0 oat 7$5x310°° 

D ore zone 15 be authnn Pe5xa07* 

ic shale 20 te1Ge7 2.0 ORF 

B shale 30 ie ei Rta 320x970" 8 

A shale 50 op a nO i 2.5x 07 > 


The storativity is assumed to be constant 
throughout at a value of 0.0001. 


For the purposes of defining the flow domain, the 


rivers bounding the three sides of the system are 
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represented by a series of constant head nodes at the 
surface. The south, bottom and top sides of the system are 
defined as no-flow boundaries. The simulation of dewatering 
the open pit mine was accomplished by defining constant head 
nodes that represent drawdown to the base of the mine. It 1s 
assumed that there is no recharge to the system and that the 
initial head distribution is consistent everywhere within 
the domain. 

To carry out this simulation, Verge's (1975) algorithm 
for modelling three-dimensional transient groundwater flow 
was used. The simulation was carried out over a period of 15 
years using a time-step of 5 years. 

The mesh generated for the simulation is illustrated in 
Figure 8. The layer shown illustrates the general 
configuration of the first of eight layers of the mesh; each 
layer having 672 nodes, 580 elements and numbered to have a 
half-bandwidth of 89. 

The mesh density was determined more on the basis of 
the memory requirements of the Verge (1975) algorithm when 
used on the Amdahl 5860 mainframe than on the Sent of 
computational efficiency or detailing the solution. 

The mesh required approximately 1.5 hours at the 
graphics terminal to input the information to construct the 
mesh. Approximately 32 seconds of CPU time, on the Amdahl 
5860 mainframe, is necessary to construct the grid 
Approximately half of this computer time was involved with 


the renumbering of the mesh to reduce the bandwidth. In 
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addition, approximately another hour of time is required to 
determine the location of constant head nodes that define 
the position of the river and the mine. 

The drawdown cone, as seen at the surface after 15 
years, is illustrated in Figure 9. This drawdown cone 
follows the outline of the mined area and provides a 
reasonable representation of what one would expect. The 
Simulation required approximately 45 seconds of CPU time on 


the Amdahl 5860 mainframe, 


von ell ai . <tugas aE 
peapeinsitng Aer eer ion exitsv0' 
oat ation Bios $ao tt vdaanoorar 


bet asad 


) 
a) 


‘HyuUeWOeO Jeu 
YGL OW! poz olosip UNPWOp OY) BUIMo ys GunBy poduys seynBuvyoos OnE ULL ounbig 


= 
; 
. 


Pe Nie 


a | . | , wot be aang 8 st 
= eet Otrsi besiiexc2ib clipe ent gniwor S1wQ “- "asnea rong stone 


NVAINVANVANANVAN 
NANNY 


ZNZ| 3 
NZNZNVNZNZN 
VANVANVANVANVANVAE 
NZNZNVNZNEN 


MANZMARVANVAN 


VANVANVANVANVANAE 
NVANVANANVANVAN | 
VANVANVANVANVANVAIE 


epomMa ise HOWSIY G~rCePSES WIOTBY 


| 


q 


ae 


ons ‘I : 
i) UM . 


ANY Nae 
ZS RARNVAISY:- iE 
AIASVANMA SS 


VARVARA 
SD ieatan. <7 
RUN — a 
ve. SNES 


pital Bs2I¢ weoraudrys spb Os apomivd ips Qowsiy giecLte1seq IO BOO MIsUANR. Sense 


ans 


: rs 
par ” 
a ° 
ha 


SKDED Kear 
aw KUN NIN ‘ 
SCIICN NS MN 


oe 7 A 
a o 


aa } 


AZEPN CNRS 
WANA ANTANVAN 2. 


sory a 
y B : 
Y we 
« 
a 


es) 
et ai 
- 

mj 

U 


4, OPTIMIZATION OF THE MESH GENERATION PROCESS 

With the availability of the mesh generation procedure, 
the problem still remains as to whether the mesh produced is 
optimum for the simulation in question. In the previous 
chapter, the results obtained from the hand-drawn mesh and 
the computer-produced mesh, for the same problem, were 
somewhat different from each other. Where it is not possible 
to verify numerical solutions with analytical solutions for 
a given problem, it is hard to know which simulation is more 
accurate. In the example considered, the computer generated 
mesh has a higher density of nodes and elements than the 
hand drawn version. However, the former could be considered 
computationally more efficient because of the smaller 
half-bandwidth. Arguably it could be stated that the 
computer produced mesh has more elements and nodes, and 
therefore should have provided a more exact solution. The 
question would then become "did we really need that many 
nodes?". 

Bathe and Sussman (1983) state that a mesh of 
relatively few elements whose nodal points are placed in 
optimal positions can match the predictive capabilities of a 
coarse mesh whose nodal points are not particularily well 
placed. As well, the use of higher-order isoparametric 
elements can improve the solution if the elements were 
distorted in such a manner that they fit the optimal 
positions of the nodes. This last point will not be 


considered further because the generation scheme presented 
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herein only generates linear elements. 

In my code two facilities were incorporated to help 
Optimize the mesh. There is an inherent advantage in dealing 
with pre-positioned flux nodes in dealing with hydrogeologic 
Simulations. In other words, because the source of the flux 
is known, it is possible to compress the mesh to maximize 
mathematical accuracy. The manner in which this compression 
is performed is illustrated in Figure 10. The existing 
elements are represented by solid lines, while the 
partitioning is represented by dashed lines. It can be seen 
that the elements surrounding a given well (or flux node) 
are partitioned to form more elements. The partitioned 
elements are then "Smoothed" using an averaging technique to 
reposition the nodes so that there is a smooth transition in 
element size away from the flux node. Time constraints did 
not allow me to investigate grid compression around a flux 
node more thoroughly. This would have involved optimizing my 
scheme and investigating other schemes such as that used by 
Bathe and Sussman (1983). 

As well there is a routine to optimize the 
computational effort involved with the solution of the 
finite element mesh. This routine renumbers the nodes in 
order to reduce bandwidth, but unlike the algorithm provided 
by Collins (1973), the routine incorporated into this 
procedure renumbers the nodes beginning at one corner of the 
generated mesh. Collins' (1973) algorithm provides the 


ability to find the optimum node numbering scheme by 
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successive renumberings, starting with a different node each 
time, and comparing the resultant bandwidths. From the 
solutions provided by Collins (1973), the best schemes begin 
numbering at a corner of the domain. It is not known whether 
this provides the optimal renumbering scheme in every case, 
but it was felt that the extra effort and cost to find the 
optimal condition would not result in savings in 
computational effort within the solution. 

Logically, the efficiency of the renumbering scheme may 
be due in part to the geometry of the mesh (ie; there is 
always a long and a short axis). It is assumed that the 
worst case Situation would be one where the mesh forms a 
large cube and the bandwidth would be related to the length 
of a diagonal through the cube. For the vast majority of 
cases, the renumbered mesh should have a smaller bandwidth 
than the original. 

Workers in the field of structural engineering have 
come up with schemes to provide optimum meshes for finite 
element simulations (Bathe and Sussman, 1983; Shephard et 
al., 1979, 1980). These procedures involve an iterative 
approach begining with a crude grid that is optimized when 
interacted with a finite element solution. In stuctural 
engineering this optimization is accomplished through 
modifying the initial mesh until the strain energy of the 
solution is maximum or the potential energy is minimized 
(Bathe and Sussman, 1983). In hydrogeology, this approach 


would be analogous to modifying the mesh until the change in 
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the potential head is minimized for each element. A 
fundamental approach to this problem is that the potential 
should be minimized not only with respect to finite element 
nodal displacements but also with respect to the nodal point 
locations (Bathe and Sussman, 1983). Defining the potential 
energy by Il, and the condition that the energy not fluctuate 


by oll=0, gives the following two equations. 


oll/dau; = 0 
ém/o6; = 0 (3) 
where: u; = the nodal point displacements 
B = all the nodal coordinates that can be changed 


without changing the geometry of the domain. 


The effect of satisfying these equations is that the 
resulting mesh is concentrated around areas of greatest 
flux. In my code there exists a facility to concentrate flux 
nodes prior to the initial simulation. Therefore, satisfying 
these equations would involve an iterative approach once the 
initial solution is performed. 

It is beyond the scope of this work to provide 
algorithms for solutions to these two equations. Indeed it 


is only recently that workers in this area have begun to 
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Figure 13. Triangular shaped figure discretized with a mesh of triangular 
elements. Corner points used to generate the mesh are 
designated with an “xX”. 
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deal with this problem (Bathe and Sussman, 1983; Shephard et 

Gat oy See hoo0 ):; 

In my work there has been a modest attempt to provide a 
degree of built-in optimization. There exist some 
approximate rules and quidelines for constructing effective 
meshes through iterative procedures that involve the 
performance of analysis. The mesh generation scheme 
presented in this work does not easily allow one to use the 
approaches outlined by such authors as Shephard et al (1980) 
and Bathe and Sussman (1983). 

It would therefore be advantageous in the synthesis of 
the mesh to try to bring in a degree of optimization before 
the simulation is attempted. Arguably this can only be 
attained through a great deal of experience, but two 
guidelines are presented with the figures to aid the 
uninitiated: 

1. Probably the most important point is that in problems of 
aquifer hydraulics, the region of greatest variation in 
hydraulic head is known. Even though the program has a 
built-in facility to handle the increased change in 
hydraulic head, the user is advised that there should be 
more elements in an area where greatest changes in 
potential occur, and less where little change is 
anticipated and they are not necessary. 

2. Although previously stated in the body of this work, 
that in dealing with isoparametric elements, it is not 


easy to change the density of the elements, there are a 
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few techniques involved in the placement of corner nodes 
that can be used to change this distribution. This is 


illustrated in Figures 13 and 14. 
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5. CONCLUSIONS AND RECOMMENDATIONS 

In accordance with the objectives set out in this 
study, a finite element mesh generator has been successfully 
developed. The mesh generation scheme has proven successful 
and expedient in modelling finite element groundwater flow 
problems presented herein. It is anticipated that the 
generator will be successful in generating meshes where the 
geometry of the domain is not extremely complicated, 
requiring only gradual changes in mesh density. It will 
accept large contrasts in hydraulic conductivity, but due to 
its inability to produce a graded mesh, the solution may be 
inaccurate. As well, it is anticipated that the scheme will 
be successful in most cases where the geometry of the domain 
boundaries do not form very acute angles (ie: generally less 
than 10 degrees). 

The amount of effort normally required to input data 
for a simulation has been reduced considerably through 
incorporating into the design of the system a digitizer and 
‘interactive graphics terminal. The digitizing tablet is used 
to specify both numerical information as well as options in 
mesh generation. The colour graphics terminal, used in 
conjunction with the Integrated Graphics facilities (Davis, 
1977) available at the University of Alberta, allows the 
user to interact in the mesh generation procedure. The 
program is capable of generating finite element meshes in 
both two- and three-dimensions and has proven to be capable 


of handling problems of the complexity commonly associated 
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with aquifer hydraulics. 

This system operates in a mode that allows the user to 
be able to generate meshes without having previously 
absorbed volumous user manuals. The program prompts the user 
with options throughout the entire generation process. This 
"user friendly" feature should allow for the wide use of 
this program. 

The graphics features display meshes in a series of 
eight colours, allowing the user to visualize where he is in 
the generation process. As well, hard copies can be obtained 
at any time during the mesh generation process. 

The program, by being uniquely developed for the 
purposes of modelling aquifer hydraulics, has the ability to 
incorporate both pumping and observation wells within the 
mesh. 

A degree of grid optimization has been incorporated 
into the program. This capability relates to two main areas, 
namely: in improving mathematical accuracy; and 
computational efficiency. To aid in the area of mathematical 
accuracy, the code compresses the mesh around flux nodes to 
better approximate rapid changes in hydraulic head around 
pumping wells. Computational efficiency is aided through the 
incorporation of the algorithm provided by Collins (1973) to 
decrease bandwidth. 

Recommended areas for further study are listed below. 
1. The area of grid optimization should be fully 


investigated. Again, as in the area of finite element 
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mesh generation, the field of structural engineering 
provides the logical starting point due to the material 
already existing on this problem in that field. It is 
suggested that the work involve using the information 
already supplied as input to the program such as 
transmissivity, storativity and pumping rates and the 
relationships of potential energy as provided by Bathe 
and Sussman(1983) as a logical starting point. 

With the advent of the personal computer and their 
increasing memory capabilities and decreasing costs, it 
would be advantageous to investigate the feasibility of 
adapting this system for use with these smaller systems 
so aS to readily allow its use by small engineering 
companies to aid in the. speedy solution of everyday 
finite element simulation problems. Problems envisioned 
in this would be adapting the program's graphics to the 
capabilities of the personal computer's ‘graphics and 


decreasing storage requirements for this program. 
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APPENDIX A User's Guide 


To use the automatic mesh generation program one should 
have at his disposal: a diagram of the figure to be 
modelled, with transmissivity zones, pumping wells, 
observation wells and areas of leakance and values of 
necessary constant head delineated; and hydrogeological 
parameters. It is also assumed that the user has the 
computer hardware or equivalent previously described. To run 
the compiled version in MTS the commands are: $RUN 
MESH.OBJ+*1IG 5=*MSOURCE* 6=*SINK* 7=-A 9=-C 
where: MESH.OBJ is the compiled version of the automatic 
mesh generation program. 

*IG is the integrated graphics library of subroutines 
available at the University of Alberta. 

5=*MSOURCE* specifies the input port (or device) that 
accepts data from the keyboard and digitizing tablet. 

6=*SINK* specifies the terminal as the output port. This 
facilitates prompting, querying and displaying the. 
results so that the user is involved with all aspects of 
the mesh generation. *SINK* displays everything on the 
terminal screen. 

=-A this is the port(or device) where the resultant mesh 
information is outputted to after the generation is 
completed. 

9=-C this is the port to which all plots requested are sent 


for plotting on the CalComp Plotter. 
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Once the program is loaded the user will be asked to 


respond to a series of queries. This section is intended to 


aid the user as to what to expect to be able to answer as 


the generation routine progresses. 


he 


IM-a variable that determines whether one is using 
English or metric units. IM=0 informs the program that 
the user 1s working in English units and the program is 
expecting to digitize a map where "1 inch=x feet". If 
IM=1, the program will convert all coordinates into 
meters. The program is expecting to digitize a map where 
the scale is defined as "1:xxxx". 

Scale-when entering the scale the format is F10.0. 
Therefore the user is cautioned to remember to input a 
period(or a comma when dealing with integer input) after 
the number is typed in such as 10. or 10,. 

ID-a variable that specifies whether one is working with 
a two-dimensional or three-dimensional mesh. ID=0 
implies that the mesh to be generated is two-dimensional 
and to be generated in the x-y plane. ID=1 implies the 
mesh to be generated is three-dimensional, generated in 
the x-y plane and projected in the z-direction. 
IELEM-prompts the user as to what type of element is 
required. In two-dimensions, the request would be for 
IELEM=0, which implies triangular elements, or IELEM=1 
for quadrilateral elements. In three-dimensions the 
choice would be between IELEM=0 for tetrahedral 


elements, or IELEM=1 for parallelpiped (or "brick-type") 
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elements. Note that this mesh generator is only capable 
of constructing linear elements. 

ISS-is a prompt only used when one is working in 3-D. 
Where ISS=0 implies Steady State and ISS=1 implies 
Transient simulation. This outputs a value of 999. for 
the time, that triggers the steady state simulation in 
Verge's (1975) 3-D Groundwater model program. 

This. prompt requests the user to indicate the "LOWER 
LEFT" "LOWER RIGHT" and "UPPER RIGHT" corners of the map 
on which the flow domain is defined. These data are 
provided so that the computer can scale the drawing for 
purposes as displaying on the screen, plotting and for 
determining whether the map is tilted with respect to 
the axis of the Pia aa ae 

The next prompt is the menu. The user is provided with 
11 options that would control the program in the manner 
outlined in the flow chart in Figure 4. The following 
prompts will be in the order as followed by this flow 
chart, and a progression through the mesh synthesis 
process. 

After selecting Switch No. 9, the user would be asked to 
enter the storativity and the transmissivity (in 3-D the 
hydraulic conductivities). The format for reading the 
values is E10.3. After this, the user would be requested 
to digitize the perimeter of the transmissivity zone in 
a counter-clockwise manner. 


If any pumping wells are present, the user would then 
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select Switch No. 7. In 2-D the prompt would be to enter 
the "Start time", "finish time", and the number of time 
steps for the simulation, and the pumping rate in 
L?/sec. The query would then be asked whether there is 
any more pumping information. This question provides the 
Opportunity to add more simulation periods as is 
compatible with the two-dimensional Verge-Frind (Verge, 
1972) transient code for confined aquifers. In 
three-dimensions the user would be requested to enter 
the simulation period,(the time units are dependent on 
the time units designated in the Verge(1975) 
three-dimensional model), initial time step, a 
multiplier for the time step, number of time steps 
changes in delta time and the pumping rate in L?’/sec. 
The format for entering the start time, finish time, 
Simulation period, initial time step, multiplier for the 
time step and pumping rate are all F10.0. The formats 
for entering the number of time steps and the number of 
time steps between changes in delta time are 14. 

The next prompt would be encountered after selecting 
Switch No. 13 for generating a mesh. The prompt would be 
to designate the four (eight in three-dimensions) corner 
points of the figure displayed on the terminal. If one 
is working in three-dimensions, a prompt would request 
the elevation of the corner point. As well there are 
prompts that request that the user double check the 


points that he has input. The next prompt would have the 
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user specify the number of elements desired in the x-, 
y- and z-directions to obtain the desired mesh density. 
Once the mesh is drawn, the user is prompted as to 
whether or not the figure displayed on the screen is to 
be plotted on the CalComp Plotter. From here,.the user 
has the option of changing the size of the plot, with 
the default being the same size as that displayed on the 
screen. 

Other Prompts include: 

Entering the leakance values in format E10.4. 


Entering. the boundary fluxes in format F10.0. 


Dimensioning of Program Arrays 


As noted in the body of the text, the arrays are all 


dimensioned sufficiently large for most situations. However 


the following is a guide should there arise a need to change 


the storage requirements. 


COMMON MEMJT(NNODE, 32) ,IN(NELM,8),IFACE(3500,6),YN(NNODE), 


ZN (NNODE) , CX(NNODE) , XN(NNODE) , JIMEM(NNODE) ,NODE1(8), 
NDX (NNODE) ,NDY (NNODE) ,NDZ(NNODE) , ITYPE(NNODE), 
NDWEL(20,2,25),LEAK(100),ALEAK(100),XIN(100), 
YIN(100),ISAVE(25,2) ,PPOSX(10),PPOSY(10),PPOSZ(10), 
OBSX(10) ,OBSY(10),0BSZ(10),ND(16),NND(16),COX(10), 
Ccoy(10),C0z2(10),IFE(6),NPMP(25),IPPN(10),ICRNR(10,8), 


NFACE(6,4),ICNTST(4,2) ,NBDND( 4000) 


DIMENSION TIME(10,2),NTIME(10),PMPRTE(10,10),X1(4),Y1(4), 
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Z1(4),TITLE(10),BFLUX(10),A(3),B(3),NCEL(10,2), 


CHNG(10),ITCHNG(10),STOR(10) 


INTEGER COL(10) 


Where: NNODE=estimated number of nodes. 


NELM =estimated number of elements. 


Use of arrays in common: 


MEMJ T 


IN 


IFACE 


XN 


YN 


ZN 


Cx 


JMEM 


NBDND 


NDX 


NDY 


NDZ 


Nec ASO 


stores the nodes connected to each node. Used 
in SETUP and OPTNUM. 
array describing incidence of each element. 
array denoting elements appearing in each of 
the six faces of the mesh. Used throughout 
the program. 
array storing x-coordinate of each node. 
array storing y-coordinate of each node. 
array storing z-coordinate of each node. 
transmissivity of each node. 
number of nodes connected to each node. Used 
in SETUP and OPTNUM. 
nodes that form the perimeter of each 
transmissivity zone. 
the adjacent node in the positive x-direction. 
the adjacent node in the positive y-direction. 
the adjacent node in the positive z-direction. 
array denoting type of each node, 

ie: O-normal 

1-constant head 


2-flux node. 
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NDWEL - defines the node numbers for observation and 
pumping wells in 2-D and 3-D. Dimensioned for 
20 wells, observation or pumping well, and 
for each well to occupy 25 nodes in a vertical 
sense. 


defines an element that a leakance value is 


LEAK 
assigned to. 

ALEAK - defines the leakance value for each element. 
XIN — f$ormanputting x-coordinates in "Pic™ routines. 
YIN - for inputting y-coordinates in "Pic" routines. 
ISAVE - defining nodes that surround a pumping well. 
PPOSX - x-coordinate of a pumping well. 

PPOSY - y-coordinate of a pumping well. 


PPOSZ - z-coordinate of a pumping well. 


OBSX - x-coordinate of an observation well. 

OBSY - y-coordinate of an observation well. 

OBSZ - z-coordinate of an observation well. 

ND - used in the concentration of elements around 
a pumping node. 

NND - used in the concentration of elements around 
a pumping node. 

COX - hydraulic conductivity in the x-direction. 

COY - hydraulic conductivity in the y-direction. 

COZ - hydraulic conductivity in the z-direction. 

IFE - used as a counter in conjunction with IFACE. 

NPMP - stores the number of nodes associated with 


each well. 
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NEACE = 


HON Si 


NODE1 - 


Stores the number of simulation periods 
associated with each well. 

stores the corner points associated with each 
transmissivity zone. 

node numbers with each face of each new 
transmissivity zone. 

store common corner points. 

book keeping array for the formation of the 


elements. 
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